Breakdown of chiral symmetry during saturation of the Tayler instability 
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We study spontaneous breakdown of chiral symmetry during the nonlinear evolution of the Tayler instability. 
We start with an initial steady state of zero helicity. Within linearized perturbation calculations, helical pertur- 
bations of this initial state have the same growth rate for either sign of helicity. Direct numerical simulations 
(DNS) of the fully nonlinear equations, however, show that an infinitesimal excess of one sign of helicity in 
the initial perturbation gives rise to a saturated helical state. We further show that this symmetry breaking can 
be described by weakly nonlinear finite-amplitude equations with undetermined coefficients which can be de- 
duced solely from symmetry consideration. By fitting solutions of the amplitude equations to data from DNS 
we further determine the coefficients of the amplitude equations. 



PACS numbers: 52.35.Py, 11.30.Qc, 07.55.Db, 47.20.Bp 

I. INTRODUCTION 

There are many examples in nature where the ground state 
does not share the same symmetries of the underlying equa- 
tions of motion |1J. The most common examples include 
equilibrium phase transition, e.g., the case of a liquid-solid 
transition where the space translational symmetry is broken, 
or that of a paramagnetic-ferromagnetic transition where the 
spin-rotational symmetry is broken; see e.g., Ref.[|2i] for a de- 
tailed discussion. However, the original symmetry is not lost 
but gives rise to the appearance of a regular structure with a 
specific length scale. 

In non-equilibrium physics spontaneous symmetry break- 
ing is often observed when some control parameter is in- 
creased above a critical value, see e.g., Ref [3] for a compre- 
hensive introduction. Two well-studied examples from fluid 
dynamics include the case of Rayleigh-Benard convection ^ 
and the Mullins-Sekerka instability of a moving interface be- 
tween two phases |4]. These systems, too, are invariant under 
translation and reflection, but the basic instability produces a 
symmetry breaking bifurcation in which the continuous trans- 
lational symmetry of the basic state is broken to a discrete 
one, although the mirror symmetry is often retained. If the 
instability parameter is raised further, secondary instabilities 
may break the periodic pattern and eventually a completely 
new symmetry-broken state may emerge, as has been seen in 
several experiments [5]. At very high values of the control 
parameter, turbulence sets in and most of the symmetries are 
statistically restored. 

In a hydrodynamic system under rotation, spontaneous 
breakdown of chiral symmetry has been studied; see, e.g., 
Ref. 101 • Spontaneous chiral symmetry breaking is also found 
in liquid crystals [6]. Preliminary evidence showing spon- 
taneous chiral symmetry breaking in magnetohydrodynamics 
(MHD), in the absence of rotation, has been presented for the 
magnetic buoyancy instability [8] and for the Tayler instabil- 
ity in a Taylor-Couette setup |9]. However the role of the 
dynamics of the bifurcation process is still poorly understood. 

The purpose of the present paper is twofold. First we 



demonstrate the occurrence of spontaneous chiral symmetry 
breaking in the context of a global instability of the toroidal 
field, and second we elucidate some aspects of the underlying 
nonlinear mechanism which determines the evolution from a 
mirror-symmetric state to a state with a preferred handedness 
or helicity. In particular, we shall be interested in the case of 
the Tayler instability iT^ITTll . which has attracted much inter- 
est in recent times for its possible astrophysical applications 
ll9l [l"2l-[T8|] . We thus discuss the possibility of generating a fi- 
nal state with finite helicity starting from a non-helical basic 
state, using a very small controlled helical perturbation. 

Our setup has the advantage of better clarifying the com- 
plex nonlinear coupling between the different modes, which 
eventually leads to the formation of a final helical state. In 
fact, the Tayler instability, in its simplest realization, has no 
threshold field, at least in ideal MHD ifTsll . where a sufficient 
condition for instability simply reads 

a In s 2 

s being the cylindrical radius. On the other hand, the spec- 
trum is characterized by an infinite number of unstable modes 
all characterized by pairs of opposite azimuthal wave number 
TO = ±1, 2, 3, but with precisely the same growth rate. In 
particular, as is well known, to = ±1 are the modes with the 
fastest growth rate. Here our aim is to understand the dynam- 
ics of the bifurcation process which leads to the selection of 
a final state of finite helicity and to understand the evolution 
of the system after the bifurcation takes place. It should also 
be noted that in the linear stage the preferred helicity is deter- 
mined essentially by the helicity of the perturbation, but the 
nonlinear evolution can be rather complex and it is not clear a 
priori what the final selected helical state would be. 

The rest of the paper is organized as follows. In Section HIl 
we write down the finite- amplitude equations that govern the 
evolution of the instability in the weakly nonlinear phase. Our 
approach is based on symmetry arguments; a detailed ana- 
lytical derivation of the amplitude equations is avoided here. 
We find that the amplitude equations predict a breakdown of 
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parity for certain choice of parameters. Direct numerical sim- 
ulations (DNS) of the fully nonlinear equations describing the 
evolution of the Tayler instability are performed in Sectionlllll 
In our DNS studies we also find breakdown of parity. We fit 
data from DNS to solutions of the amplitude equations to nu- 
merically determine the parameters appearing in the ampli- 
tude equations. It turns out that the amplitude equations we 
deduce are identical to those used to describe the breakdown 
of mirror symmetry in studies of the biochemical origin of life. 
This connection is explored in Section|IV] Finally conclusions 
are drawn in Section IV] 



II. AMPLITUDE EQUATIONS 



where the simplest form of the Lagrangian is given by 



C[L,R] =7 



\R\- 



(8) 



The form of the Lagrangian is determined by the symmetry of 
the problem. The coefficient 7 is the linear growth rate and ji 
determines the saturation of the instability in the weakly non- 
linear regime. We emphasize that the /i and 7 for L and R 
could be different if and only if the chiral symmetry is broken 
from the outset, but this is not the case here. Now note that 
the Lagrangian must also be invariant under the parity trans- 
formation, under which 



V{L) = R and V{R) = L. 



(9) 



To the best of our knowledge the amplitude equations de- 
scribing the spontaneous breakdown of mirror symmetry in 
hydrodynamic instabilities were first described in Ref. I119ll . 
The basic idea is as follows. 

Let us consider an instability with two growing modes with 
opposite helicity but exactly the same growth rate and let the 
amplitude in this basis of the left- and right-handed modes be 
given by vectors L and -R, respectively. In physical space we 
have 



L(a;) = L(j){n), 
R{x) = R(j){n). 



(2) 

(3) 



For example, in Cartesian domains, with real-space coordi- 
nate X, 0(n) = exp(m • x). In cylindrical coordinate, is 
a combination of trigonometric and Bessel functions. As the 
modes are helical they satisfy the Beltrami relation. 



V X R = AR and V x L = -AL. 



(4) 



For 



the present problem, explicit expres- 
sions involve a linear combinations of the type 
Jm(,s^J A? + n^vr^/Zi^) cos(m^) cos(zn7r//i), where Jm 
is the Bessel function of the first kind, n,m = ±1, 2, 3 . . ., h 



is the height of the cylinder and s the cylindrical radius |12C 
The set of such modes forms a complete set (a Hilbert basis) 
for the spatial distribution of the field. 

Here we assume that the dynamical evolution of the unsta- 
ble mode is determined by an effective Lagrangian. For the 
left-handed helical mode, total helicity and energy are given 
by 



El = \j h\x)<fix = H-L* 



-Hi 



-2AEr 



(5) 



(6) 



where an asterisk denotes complex conjugation. Analogous 
definitions apply also to Eji and Hr — +2AEii. We then 
have E = El+Er being the total energy and H = Hl+Hr 
the total helicity. In the weakly nonlinear regime the ampli- 
tude equations can be written as 



dL SC 
— — — and 
dt SL 



dR_5C_ 
dt " 5R' 



i1) 



This additional symmetry allows one additional term in the 
Lagrangian given by 



L?\R[ 



(10) 



With this additional term in the Lagrangian the evolution 
equations for the two eigenmodes are given by 



9i 

'dt 
OR 

'dt 



= 7i- (M|iP+M*l-Rp)i, (11a) 

==7^- (mI-RP + m*!^!')-^- (lib) 



These equations, for certain parameters, allow and can de- 
scribe the growth of one handedness while the other is ex- 
tinguished [19(]. Similar equations, which describe the time 
dependence of the amplitudes of the leading modes, but with- 
out considering their spatial dependence, are often used to ex- 
tend linear perturbation theory of hydrodynamic instabilities 
into the weakly nonlinear regime. In this form they are often 
called the Landau equations 1I2III . 

The energy of the left and right-handed modes is deter- 
mined by, 

''^^ : 2-fEL-A^iEl-A^i,ELER, (12a) 



dt 

dER 

dt 



2jEr - A^iEl - ^pl.ElEr. 



(12b) 



These equations show that both E^ and Er_ grow exponen- 
tially at the rate 27 until nonlinear effects become important 
and either El or Er saturates at i?o = 7/2/i and the energy 
of the mode of opposite handedness vanishes. In principle, the 
achiral solution with E^ — Er = E'a = 7/2(/i + /i^) is also 
possible, but, as we will see below, such a solution is unsta- 
ble for ^ < fi^,, which is what we find Sect. |III1 The reason 
for this instability is the presence of the term proportional to 
/i*, which represents a phenomenon known as "mutual an- 
tagonism" in studies of the origin of homochirality of bio- 
molecules ll22l - l24ll . We will return to this issue in Sect. |Iy] 
where we discuss the analogy with chiral symmetry breaking 
in biomolecules in more detail. 

Using Eqs. ^ and ( fTTT i and defining H — 'H/2A we have 
Er = {E + H) /2 and El^{E- H) /2. We can thus obtain 
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the following evolution equations 



dE 



^ = 2-fE-2{fi + n,)E^^2{fi-fi,)H^, (13a) 



^ = 2-/H-AnEH. 
at 



(13b) 



The dynamical system described by (fT3T l and depicted in Fig.[T| 
has four fixed points in the {E,H) plane. Si = (0,0), 
52.3 = {Eq,±Eo), and Si = (2Ea,0) with eigenvalues 
Ai = (27,27), A2 = A3 = (-27,2(/.t - /i,)/7), and 
A4 = (—27, 27 — 47/i/ (/i + /i*)). The origin is always re- 
pulsive while 52 and ^3 are sinks or saddle points depending 
on the values of parameters /i and /i,. ^4, corresponding to 
the achiral solution, can be an attractive point only if /i, < /i, 
otherwise is a saddle point. 

A discussion of the amplitude equations is now in order 
Firstly, we have assumed that there are exactly two modes of 
opposite helicity that becomes critical at the onset of the insta- 
bility. This assumption is based on linear perturbation analy- 
sis. As all the other modes in this case stable, in the spirit 
of center manifold reduction, we have ignored their contribu- 
tions to total energy and helicity. If several modes are simul- 
taneously unstable at the onset, then we may expect a higher 
degree of complexity. Secondly, as our approach is based on 
symmetry, the form of the amplitude equations that we obtain 
is very general. This is also a weakness of our approach, as 
we cannot determine the expression for either fioi fi^,. In prin- 
ciple, the method of multiscale expansion or center manifold 
reduction can be applied to this problem to derive analytical 
expressions of /i and /i*, but this is a difficult proposition in 
the present case as a solution of the linear problem itself is not 
known in an analytically closed form. 

We can then compute the quantities 7, /i and /i, with the 
help of direct numerical simulations (DNS) by comparing the 
time evolution obtained for the left-hand side of (l3[ in the 
weakly nonlinear phase where our description is valid. We 
can anticipate that in most of our simulations /i < /i* and 
therefore the system should relax to a state of finite helicity in 
a finite time, although we start from an infinitesimally small 
hehcity. This is precisely what we observe in our DNS. 



III. DIRECT NUMERICAL SIMULATIONS 

To analyze the evolution of the Tayler instability we choose 
our numerical domain to be a cylindrical shell with an inner 
radius Sin = 1, outer radius Sout = 3, and height h = 2. We 
perform simulations of the time-dependent resistive magne- 
tohydrodynamic equations for a compressible isothermal gas: 
the pressure is thus given hy p ^ pc^, where p is the density 
and Cs is the isothermal sound speed. 

We use the PENCIL Code [SOJ to solve the equations for 
the magnetic vector potential A, (B = V x A) the velocity 
U, and the logarithmic density In p in the form 



dA 

'dt 




FIG. 1: The phase portrait for /i < /i,. This is the typical situation 
in which 52 and 5*3 are attractive and 5*4 is a saddle point. 



^ = -c^ V In p + J X B/p + 



Dln/9 
Bt 



-V-C/, 



(15) 



(16) 



where 



-Fvisc = P ^ V • 2vpS 



is the viscous force, S is the traceless rate of strain tensor hav- 
ing components Si^ = \{Ui^j + Uj_i) — \5ijW ■ U, 

J = V X B/po 

is the current density, ly is the kinematic viscosity, and i] is the 
magnetic diffusivity. 

We choose periodic boundary conditions in the vertical (z) 
and azimuthal (tp) directions, while at radial (s) boundaries we 
select perfectly conducting boundary condition for the mag- 
netic field and stress-free boundary conditions for velocity. 
The resolution of the simulations presented here is 128"^ mesh- 
points in all three directions, but comparison with different 
resolution demonstrated that our results are converged. 

We choose a basic state with zero velocity and zero axial 
component of the magnetic field (Bz). The azimuthal compo- 
nent of the magnetic field is 



= Bo (s/so) cxp[-(s - so)Vct^ 



(17) 



U xB- 



(14) 



where Bq is a normalization constant, sq = 2 and cr ~ 0.2. 
We choose Bq and Cg in such a way that the sound speed 
is much larger than the Alfven speed. In this way we avoid 
magnetic perturbations to be dominant over hydrodynamical 
perturbations. 

In the basic state the Lorentz force due to the magnetic field 
is balanced by the gradient of pressure. Hence the pressure of 
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the fluid is given by 



„2 ,2^„-2 ^"-'f 



{2s'-cj')e- 



(18) 



where is a constant that must be large enough to ensure 
that the pressure is positive. If no perturbation is added, the 
system remains stationary. Therefore, we add at the beginning 
of the simulation a perturbation of the magnetic field with an 
infinitesimally small net helicity given by the following ex- 
pression: 



A = 5s cos 



sm nriLp 


.cos mif j 



(19) 



where 5 is an arbitrary small amplitudes which we set to 10~^ 
for all the simulations and fc^ ~ q/ = nzir/h is the vertical 
wavenumber of the perturbation. 

As discussed in |13] kink instabilities are special case of 
the so called quasi-interchange instabilities, where combined 
azimuthal and vertical field are present in the basic state. In 
the incompressible limit the unstable eigenmodes can be de- 
scribed by a (t, z, (p)-dependence of the type oc cxpijt — 
ikzZ — irrup) where the growth rate 7 is determined from a 
numerical solution of the nonlinear eigenvalue problem for 
the radial disturbance vis 



d 
ds 



1 / 2 , dvu 



m(l + A) 



I'l. 

+ — 

s 

1 + A 

Wis - 



(cjAz + 2mu}B) (20) 



0. 



Here oja — {B ■ k)/y/p with k — {0,m/ s, kz), so ujaz = 
kzBzl^. Furthermore, we have defined ujb = B^p/s^ 
and and A = 1 + m? /s^kl. 

Equation ( l20b describes the stability problem as a non- 
linear eigenvalue problem. This equation was first derived 
by Freidberg ll25ll in his study of MHD stability of a dif- 
fuse screw pinch (see also flS]). The author found that, for 
a given value of kz, it is possible to obtain multiple values 
of the eigenvalue 7, each one corresponding to a different 
eigenfunction, and calculated 7 for the fastest growing fun- 
damental mode. The most general form of Eq. (l2Qt . taking 
into account compressibility of plasma, was derived by Goed- 
bloed flS]. Since we study the stability assuming that the 
magnetic energy is smaller than the thermal one, the incom- 
pressible form of Eq. (|20] | can be a sufficiently accurate ap- 
proximation, as we have verified. In our the case at hand, 
[B ■ k)/ y/p = mujB = because we are interested in pure 
kink (Tayler) instabilities, with Bz = in the basic state. Note 
that as uJAz = in our case, Eq.(l20li is invariant for m — )■ — m. 



In this latter case, once (|20] | is solved and vis is obtained, 
the expressions for the other perturbed quantities denoted by 
the subscript "1" , read 



Bis — B^p Vis, 

7s 



B^ 



Blip — Bip viip (/3 - 

7s 7s 



l)vis 



Biz 



Vl^ 



7s 



-B^ Viz 



d 



{kzs)^X ds 



(svis) 



2imujg Vis 



id, . m 



(21a) 



(21b) 



(21c) 



(21d) 



(21e) 



Unfortunately even for the case of pure kink instabilities, ( l20b 
cannot be solved analytically and one has to determine the dis- 
persion relation numerically. Therefore, to test our numerical 
setup we have solved numerically Eq. (|20] | for the basic state 
( [TtI i for various values of Bq and a in the limit of small va/cs 
ratio to check that in the linear phase the growth rate extracted 
from the DNS is in agreement with the linear theory. In par- 
ticular, as the inner radius of the cylinder is not at s = 0, we 
have set vis = at both inner and outer boundaries. Note 
here that the growth rate and eigenfunctions of this instabil- 
ity are known for the ideal MHD limit. Hence to compare 
with those results, we choose viscosity and magnetic diffusiv- 
ity such that the dissipative time scales are much larger than 
the characteristic growth time (inverse of 7) of the instability. 

The results are shown in Fig. |2]for the dimensionless growth 
rate V = ^tA, where Ia = Souty/p/ Bq is the Alfven travel 
time, as a function of the dimensionless vertical wave number 
q = kzSin for model Held in Table U with Bq = 0.5, and 
Uz — 10. In particular to compare the growth rate obtained 
from our DNS we have determined the characteristic verti- 
cal wave number of the unstable mode in the linear phase by 
means of the Fourier analysis of the magnetic fields. We also 
found that in all the simulations the azimuthal wave number 
of the fastest growing mode turned out to be always m — ±1 
as higher values of |m| have a smaller growth rate as shown 
in Fig. |2]l. We found that the corresponding growth rate de- 
termined from the linear phase of our direct numerical simu- 
lation is about 7 — 5% smaller than the linear value, we think 
this acceptable in view of unavoidable numerical diffusion in 
three-dimensional numerical simulations. 

We see that the eigenfunction is rather localized for q ^ 1, 
as visible in the example shown in Fig. [3] We can exploit 
this property to obtain approximate explicit expressions for 
the growth rate at large values of q. In fact we can consider 
the magnetic field approximately constant around s — sq and 
apply the small-gap approximation (see ifisll for details) so 
that Vis oc sin(7r(s — so)/f) and the dimensionless growth 
rate reads 



= - 



+ 



2cA^m^ {(J3 - l)m2 + {(3 - 3)q'^) 
2(/3 - l)A2g2 _ c2m2 (A2 (m^ - 3(?2 

A2 (m2 + g2) 4. 7^2 



(22) 
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FIG. 2; The dispersion relation for the dimensionless growth rate 
r for the m — itl mode (solid line) and for the m — ±2 mode 
(dashed). Higher values of \m\ have even smaller growth rates. This 
curve is obtained for a linear model with physical parameters corre- 
sponding to the nonlinear model Held, for which we indicate, with 
a rhombus, the growth rate for its faster growing mode. 



where c = B^/Bq ss const and A = 2a j sq. In the limit 
\, despite the uncertain approximation that we have per- 
formed, expression (l23T l differs by only 20% from the numer- 
ical solution. 

It is interesting to notice that, by using (ISTT l in the limit 
g ^ 1 we obtain the explicit expressions 



• V X ; 
(Bi-VxBi) 



4mi?^..|(^;L) 



(23) 



(24) 



where the symbol (•) denotes volume averaging. It is therefore 
clear that, at the nonlinear level, eigenfunctions with non-zero 
TO will produce both kinetic and magnetic helicity whose sign 
will depend on the sign of to. The relevant point is that modes 
with opposite m have identical growth rate, but opposite ki- 
netic and magnetic helicity and the fate of the final helicity is 
decided by the competition of modes with opposite azimuthal 
wavenumber. 

Moreover, according to the oscillation theorem frf\ as the 
771 = ±1 are unstable, all the other modes with m = ±a, 
where a > 1 is a positive integer, are also unstable, but with a 
smaller growth rate. As a consequence, although in the linear 
phase the m = ±1 modes dominate the linear growth, already 
in the weakly nonlinear phase the contribution of modes with 
TO 7^ ±1 can be also important for the selection of the final 
helical state. 

The eigenfunctions appear clearly in our simulation and 
they fit quite well the eigenfunctions calculated by the lin- 
ear model, as shown in Fig. [3] In our simulations, during the 
growing phase of the instability we observe a net increase of 
the helicity, as shown in Fig. |4] where we plot the time series 
of the normalized kinetic, current and magnetic helicity. It is 
interesting to notice that while kinetic helicity decays on the 




FIG. 3: Eigenfunction v\s for the m = 1 mode for q — 16. The 
result of the simulation, model Held in Table |I] (solid line) is over- 
plotted on the eigenfunction obtained solving ( I20t (dashed line); see 
I13n for more details. This is observed at t/tA = 9, that is during the 
linear growth of the instability. 



viscous time scale, the current and magnetic helicities reach a 
nonzero value at very large times. 

In actual simulations we choose ly = lO^'^ (in code units), 
so that the viscous time scale is t^, — jv ^ and the 
actual value of v does not play a significant role in the weakly 
nonlinear phase as we verified in our simulations. Moreover 
we decided to use a very small value for the magnetic dif- 
fusivity, 77 — 10~^ (in code units). This is done to prevent 
the decay of the magnetic field by diffusion. In general such 
small values of magnetic diffusivity would imply extremely 
large values of magnetic Reynolds number which would be 
impossible to resolve with the resolutions we use. Neverthe- 
less we choose such values to have a toroidal field stable on 
time scales much longer than those of the instability. How- 
ever, in our simulations no sharp gradients of the magnetic 
field develop, which is the reason why such small values of 
magnetic diffusivity are permissible. 

We can now determine the coefficients 7, [l and /i* using 
the time evolution of Bit) and E(t) obtained with our DNS 
in solving the model (fTsl l. This can be done via a direct two- 
parameter minimization because the exponent 7 can easily 
be determined from the linear evolution and one is left with 
only the determination of and /i,. An example of this ap- 
proach is depicted in Fig. |5] where the agreement with our 
numerical simulations is explicitly shown. Note that around 
tjtA ~ 6.5 we enter the deep nonlinear phase and our treat- 
ment does not apply anymore. We estimate this cutoff time 
for our simulations to be in the middle of the decay transition 
for d\nB/dt and d\nE/dt depicted in Fig. |5] and we have 
checked that the values of fi and /i, are not strongly depen- 
dent on our cutoff time. 

Our results are summarized in Table |I] We see that the coef- 
ficients /i and /i, are unchanged for models that differ only in 
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FIG. 4: (Color online) Kinetic, current, and magnetic helicities for 
two different runs (Models Hel and He 1ml in TableUl with helical 
perturbation and m = ±1. t is in units of the Alfven travel time Ia- 
The viscous time is ~ 10^ * Ia and the magnetic diffusion is ~ 
10® * tA- It is clear the difference in the evolutions of the kinetic, 
current and magnetic helicities, in the first, second and third panel. 
These plots show how these quantities grow with the same rate, but 
different sign, depending on the sign of the initial perturbation, that 
is the sign of m. Note that for each model the magnetic helicity has 
opposite sign of the kinetic and current helicities. 

the sign of m in the perturbation. This is what we expect and 
one of the symmetries we have used to write the Lagrangian 
(O. Model Helc shows that the growth rate depends on the 
value of Cg, but this does not change the values of and /i*. 
Model H e 1 b and model Held have a smaller growth rate due 
to a smaller va/cs- Helb has /i and /i* smaller than Held 
, due to the fact that in the latter modes with higher have 
been excited by the initial perturbation. Note that in our setup 
the ratio va/cs depends on Bq, but not on c^. This is due 
to the fact that the model is isothermal and the initial radial 
balance is obtained through a pressure gradient that balances 
the Lorentz force. An increase of riz of the perturbation, as in 
model HelnlO, leads to a similar growth rate, but smaller fi 
and This can be explained saying that, while in the linear 
phase this model evolves similarly to any = 1 model, in 
the weakly nonlinear phase the evolution is different because 
of a faster growth of modes with higher kz- In our models we 
measure 2.2 < < 2.6 . 

IV. HOMOCHIRALITY IN BIOMOLECULES 

It is instructive to consider Eqs. (fT2T i as evolution equations 
for the concentration of two molecules of opposite handed- 
ness, L and R. Let us assume that L and R can be synthesized 
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FIG. 5: Time evolution for the logarithmic derivative of kinetic en- 
ergy (solid line) E and kinetic helicity H (dashed) as measured in 
DNS for models Hel and He 1ml (see|ll(. t is in units of the Alfven 
travel time tA- We overplot a fit of the model with equations 1131 
The best fit is obtained for 7 = 2.71/tA, /i = 7.5 • tyi/sout and 
fit — 18 ■ tA/Saut and the solutions are over-plotted on the DNS 
results. 



TABLE I: For every model s goes from 1 to 3, z from -1 to 1, the 
perturbation has an amplitude S — , a — 0.2. 



Model 




Va/Cs 


Cs m 




1-tA ^^4f 


tA 




Hel 


10-1 


0.3 


10 -1 


1 


2.71 


7.5 


18 


2.4 


Helml 


10-1 


0.3 


10 +1 


1 


2.71 


7.5 


18 


2.4 


Helc 


10-1 


0.3 


20 -1 


1 


6.2 


7 


18.5 


2.6 


Helb 


2.5 ■ 10-^ 


0.15 


10 -1 


1 


2.2 


1 


2.3 


2.3 


Held 


2.5 ■ 10-^ 


0.15 


10 -1 


10 


2.2 


3 


7.3 


2.4 


HelnlO 


10-1 


0.3 


10 -1 


10 


2.75 


4.5 


10 


2.2 



from a substrate S through auto-catalytic reactions of the form 

S L, S A R. (25) 

Autocatalytic reactions of this type have been confirmed in 
laboratory experiments lEsll . Let us furthermore assume that L 
and R are capable of polymerizing to form homochiral dimers, 

L + LALL, R + RARR, (26) 

as well as heterochiral dimers, 

L + R LR, (27) 
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then the evolution equations for the various concentrations are 



rf[S] 
dt 

m 

dt 

dt 
d[LL] 

dt 
d[RR] 

dt 
d[LR] 

dt 



which obeys the conservation law |24 

[S] + [L] + [R] + 2[LL] + 2[RR] + 2[LR] const. (29) 

These equations represent a subset of a more general poly- 
merization model 1I23I1 . Comparing with Sect. |II]we see that 
Eqs. ( |28b| ) and ( |28c| i are identical with Eqs. jllaj and ( |12bt 
when substituting [L] — El and [R] = En, and identifying 



7 rci / n 1 


1 TDi ^ 

+ [RJ): 




(28a) 


7 rci n 1 


07 n 1 2 
2ks[\-\ - 


- 2fc/[L][R], 


(28b) 






-2A:,[L][R], 


(28c) 


ksn\ 






(28d) 


ksn\ 






(28e) 








(28f) 



kc = 27, ks = 2^, fc/ = 2^* 



(30) 



Hence we demonstrate quantitatively that the spontaneous 
production of helicity from the fully nonlinear system of hy- 
dromagnetic equations can be described by the simple model 
equations ( fT2] i. which in turn represent a simple set of chemi- 
cal reactions (l25Tl-(l27Ti. 

The analogy with homochirality in biochemistry is useful 
because it helps identifying the phenomenon of mutual antag- 
onism as the main cause of chiral symmetry breaking. This 
effect corresponds to a contribution to the nonlinear terms 
that result from the interaction between modes of opposite 
handedness. These are the terms proportional to /i, and fc/ 
in Eqs. ( fT2] i and ( l28T l. respectively. In the synthesis of polynu- 
cleotides this is known as enantiomeric cross-inhibition and 
has been identified in laboratory experiments ll29l . The syn- 
thesis of heterochiral dimers is essential in that it corresponds 



to the production of waste needed to eliminate building blocks 
of that handedness that is later is to disappear completely. 



V. CONCLUSIONS 

We have shown how net helicity is produced by the addi- 
tion of a small helical perturbation to a non-helical system, 
thus driving the system to a final state characterized by a finite 
value of the helicities and, therefore, breaking the initial sym- 
metry. We have shown further that this spontaneous symme- 
try breaking can be described by weakly nonlinear amplitude 
equations ( fT3] l. Furthermore, we have numerically determined 
the coefficients appearing in the weakly nonlinear amplitude 
equations ( flj] ) for the Tayler instability. Direct numerical sim- 
ulations show that the ratio between the coefficients describ- 
ing the weakly nonlinear phase is almost constant. The agree- 
ment between the analytical model and the numerical solu- 
tions is rather good in the beginning of the weakly nonlinear 
phase, as shown in Fig. |5] This demonstrates quantitatively 
the close analogy between helicity production in hydromag- 
netic flows and the development of homochirality in biochem- 
istry, which is described by the same system of equations as 
those resulting from the amplitude equations of the weakly 
nonlinear model of the Tayler instability. It will be useful to 
extend our analysis by means of a Landau-Ginzburg descrip- 
tion of the amplitude equation by including non-homogeneous 
term in our Lagrangian to discuss the possible pattern forma- 
tion in this type of spontaneous chiral symmetry breaking. We 
hope to address this issue in a following communication. 
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